Nonholonomic systems are mechanical systems subject to unintegrable differential constraints in the sense of Forbenius theorem, which have been attracting a great deal of attention in mechanics, physics, and engineering.^[1–3] Recently, the study of nonholonomic mechanics focuses on two aspects: one is European scholars who developed the theoretical framework of geometrodynamics for nonholonomic systems,^[4–8] and the other is North American scholars whose achievements mainly embody the symmetry reduction of the nonholonomic system and its application in motion planning and control of nonlinearity.^[9–12] Modern differential geometry, such as fibre bundle, symplectic manifold, and Poisson manifold, has been fully used in many aspects.^[13,14] People turned their attention to the geometric numerical integration for nonholonomic systems driven by the researches of reachability and controllability of nonholonomic systems, motion and route planning, optimum control, path tracking, stabilization of dots and sets.^[15–23] With the development of theoretical verification and application of nonholonomic systems, geometric numerical integration is not only a powerful supplement for geometrodynamics of nonholonomic dynamics, but also a connection between the theoretical research and practical application.^[24] However the works of geometric numerical integration of nonholonomic systems are limited^[18–26] and mainly focused on extension and transformation of the existing structure-preserving algorithms, such as symplectic-preserving algorithm and Lie group algorithm. The discrete variational methods have been adopted to discrete directly the d’Alembert–Lagrange principle, and obtain the discrete Euler–Lagrange equations that can preserve the momentum or energy conservation.^[18–22,26] Moreover, these methods are only suitable for nonholonomic systems subject to linear constraints, and the energy is still a conserved quantity for the system that can be transformed into Hamilton or Lagrange systems, but it is useless for the system subject to the general nonholonomic constraints.

In nonholonomic Lagrangian or Hamiltonian mechanics, the constraint submanifold of system is not a symplectic manifold, and the symplectic form constructed from the Lagrangian or Hamiltonian is no longer preserved as in the unconstrained case. Moreover, in the case of a nonholonomic system with symmetry, the momentum map is not conserved in general, due to the presence of the constraint force.^[19] So the geometric numerical integrator for nonholonomic systems is not a symplectic-preserving algorithm. For some special nonholonomic system satisfying the self-adjoint conditions, we can transform the integral problems of nonholonomic motion equations into integral problems of conditional holonomic equations. The symplectic algorithm can be used to solve the special nonholonomic problems.^[24,26] But for most of general nonholonomic systems, the Lagrangian or Hamiltonian formulation for a dynamical system is not directly universal if the physical variables remain without using Darboux transformations, and the equations of motion are not self-adjointness. Based on the Cauchy–Kowalevski theorem of integrability conditions for partial differential equation, we can reduce the equations of motion of nonholonomic systems to a first-order form, which satisfy the conditions of self-adjointness, such as the Helmholtz’s conditions.^[27] Because of the existing nonholonomic constraints, the first-order form is a generalized Birkhoffian system whose regular conditions are broken.^[27,28] Then we can construct the structure-preserving algorithms for universal nonholonomic systems in the generalized Birkhoffian frameworks. The systems whose equations of motion are represented by the first-order Birkhoffian equations on a presymplectic or a contact manifold spanned by the physical variables are called generalized Birkhoffian systems, and they are self-adjointness. So the new algorithm for Birkhoffian systems can be named as a self-adjoint-preserving algorithm.

In Section 2, the generalized Birkhoffian system and its geometric properties are introduced, and the generalized Birkhoffian realization of nonholonomic systems is given. Section 3 constructs the discrete variational algorithm for generalized Birkhoffian systems and gives out the generalized discrete Birkhoffian integrator. In Section 4, an example is given to demonstrate the effectiveness of the numerical algorithm proposed. Conclusions are drawn in Section 5. For convenience, if there is no special illustration, Einstein’s summation convention is used throughout this paper and α, β = 1,2,…,l;i,j,k = 1,2,…,n; I,J = 1,2,…; m = 2k + l; μ,ν = 1,2,…,n − l.

Consider the generalized Birkhoffian system on the m-dimension smooth manifold M with local coordinates {t,a}. Because the value of m is not always even, the manifold M is a pre-symplectic manifold, and the configuration space of the generalized Birkhoffian system is a pre-symplectic space. Denoting the generalized Birkhoffian B : M → R and functions R_I : M → R, the generalized Pfaffian is

3.1. Generalized Discrete Birkhoffian equation and its geometric properties

In order to study numerically the generalized Birkhoffian system, the pre-symplectic space is discretized firstly. With discrete coordinate , define the discrete configuration space as M × M and the generalized discrete Pfaffian P_d : M × M → R as

where k is the discrete points, t_k is the discrete time coordinates, and a_k is the discrete space coordinates. Take c_d ∈ C_d as the discrete curves in discrete space. The discrete curves are sets (t_k,a_k) ∈ M of discrete points (t_k,a_k) ∈ M. Define the discrete path space as

Definition 1 Generalized discrete action functional is a map A_d : C_d → R with coordinates

Then we present the generalized discrete Pfaff–Birkhoff principle.

Theorem 1 There exist a unique C^k−1 (k > 1) mapping D_B : M × M → T*M and unique one-forms and for the C^k generalized discrete Pfaffian P_d : M × M → R in the discrete state space M × M. For the discrete variation of all curves δc_d ∈ TC_d, we have

where map D_B is called the generalized discrete Birkhoffian map and one-forms and are the generalized discrete Birkhoffian one-forms, which have coordinate expression

The theorem can be proved by using the contemporaneous variations methods.

It is clear from Eq. (18) that c_d is a solution if and only if dA_d (c_d) · δc_d = 0 at all points expect the endpoints 0 and N. This statement at k reads

Equation (22) is named as the generalized discrete Brikhoff equation, which can be separated into the time and configuration components

Because the generalized discrete action is a functional of the entire paths, not just some discrete associated curves, the discrete paths satisfying Eq. (23) do not automatically satisfy Eq. (24). That is, both Eqs. (23) and (24) are necessary in solving the generalized Birkhoffian equations numerically and contribute restrictions on the space of solutions.

Equation (22) defines a generalized discrete Birkhoffian map F_Pd : M × M → M × M in coordinate

If the points satisfy map (25), the discrete curve c_d ∈ C_d is a solution of the generalized discrete Birkhoffian equation (22) for each k = 1,···,N − 1. There are two one-forms (20) and (21) which arise from the boundary terms in the discrete case. By using d² = 0, we have

Defining generalized discrete action functional (17) as a map A_d : M × M → R in coordinates A_d (t₀,t₁,a₀,a₁) = A_d (c_d) with the generalized discrete Pfaffian (16), the discrete solution space C_Pd can be taken as the initial condition space of system M × M. Here c_d ∈ C_Pd is the solution of generalized discrete Birkhoffian equations with the initial conditions c(0) = (t₀,a₀) and c(1) = (t₁,a₁). From the variation (18), we can prove that the one-forms (20) and (21) satisfy

With the fact that d²A_Pd = 0, we obtain the conservation law

where are the generalized discrete Birkhoffian two-forms as well as the discrete Birkhoffian pre-symplectic forms. Take an arbitrary subinterval of the discrete sequence 0,1,···,N and perform the same procedure, we can obtain the conservation law in the discrete interval. So the generalized discrete Birkhoffian map is a generalized discrete pre-symplectic map.

Consider a system subjected to a nonholonomic constraint

If we take a¹(0) = 1, a²(0) = 1, a³(0) = 0 as the initial conditions, Eq. (33) has analytical solution

Figure 1(a) gives out the analytic motion curve in the plane of (a¹,a³).

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The motion curve of system subjected to q̇1 + tq̇2 − q2 + t = 0, obtained by analytic solution (a), Bikhoff’s framework (b), and Runge–Kutta method (c).

Now we adopt the discrete variational method to solve numerically Eq. (39) in the generalized Birkhoffian framework. According to Eq. (15), the discrete Pfaffian can be obtained by using the Euler mid-point scheme^[10]

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Errors of (a) q2(a1) and (b) q1(a3) computed by discrete variation method (solid line) and Runge–Kutta method (dotted line).

In this paper, we study the numerical method of the general nonholonomic system in the generalized Birkhoffian framework. When the motion equations of nonholonomic system are transformed into the generalized Birkhoffian equations, the new equations satisfy the conditions of self-adjointness. By using the discrete variational method, we construct a numerical method for the generalized Birkhoffian equations called the self-adjoint-preserving algorithm. With the numerical example, we conclude that the discrete variational method in the generalized Birkhoffian framework can preserve perfectly the inherent quality of the original dynamical system, and relatively accurate numerical results have been obtained.